#front matter
rm(list=ls())
source('CARproper.R')

###LOAD DATA###
ewm<-read.csv("ewm.csv")


###CREATE WEIGHT MATRIX###
N<-46
num=c(4, 4, 6, 2, 6, 3, 3, 2, 5, 
5, 5, 4, 5, 3, 7, 3, 1, 4, 
5, 3, 4, 4, 7, 4, 3, 3, 
5, 5, 4, 3, 5, 6, 3, 6, 2, 2, 
5, 8, 4, 5, 3, 5, 2, 5, 4, 5)
adj = c(
38, 22, 9, 8, 
40, 27, 5, 4, 
39, 38, 32, 23, 22, 16, 
33, 2, 
46, 40, 32, 27, 14, 2, 
35, 28, 19, 
34, 26, 18, 
9, 1, 
38, 36, 29, 8, 1, 
46, 43, 40, 33, 24, 
45, 23, 15, 13, 12, 
31, 20, 15, 11, 
45, 37, 23, 21, 11, 
32, 23, 5, 
44, 42, 38, 31, 23, 12, 11, 
39, 22, 3, 
25, 
44, 42, 34, 7, 
41, 35, 28, 25, 6, 
45, 31, 12, 
45, 37, 30, 13, 
38, 16, 3, 1, 
38, 32, 15, 14, 13, 11, 3, 
46, 37, 30, 10, 
41, 19, 17, 
34, 28, 7, 
40, 39, 32, 5, 2, 
41, 34, 26, 19, 6, 
42, 38, 36, 9, 
37, 24, 21, 
44, 34, 20, 15, 12, 
39, 27, 23, 14, 5, 3, 
43, 10, 4, 
44, 31, 28, 26, 18, 7, 
19, 6, 
29, 9, 
46, 30, 24, 21, 13, 
42, 29, 23, 22, 15, 9, 3, 1, 
32, 27, 16, 3, 
46, 27, 10, 5, 2, 
28, 25, 19, 
44, 38, 29, 18, 15, 
33, 10, 
42, 34, 31, 18, 15, 
21, 20, 13, 11, 
40, 37, 24, 10, 5)
B<-matrix(0,nrow=N,ncol=N) 
k<-0
for(i in 1:N){
	for(j in 1:num[i]){
		k<-k+1		
		B[i,adj[k]]<-1	
	}
	}
ewm.W<-B/apply(B,1,sum)


#run models
ewm.1<-CAR.Proper(X=cbind(ewm$ideol), y=ewm$aveelite, W=ewm.W, mcmc=100000, beta.mean=c(.76),beta.var=diag(c(1/100),nrow=1,ncol=1))
ewm.2<-CAR.Proper(X=cbind(ewm$ideol,ewm$aveelite), y=ewm$pid, W=ewm.W, mcmc=100000, beta.mean=c(.83,-1.01),beta.var=diag(c(1/30.864,1/30.864)))
ewm.3<-CAR.Proper(X=cbind(ewm$ideol,ewm$aveelite,ewm$pid), y=ewm$demstr, W=ewm.W, mcmc=100000, beta.mean=c(.15,-.24,.78),beta.var=diag(c(1/51.020,1/44.44,1/100)))
ewm.4<-CAR.Proper(X=cbind(ewm$aveelite,ewm$demstr), y=ewm$legideol, W=ewm.W, mcmc=100000, beta.mean=c(1.05,.63),beta.var=diag(c(1/625,1/625)))
ewm.5<-CAR.Proper(X=cbind(ewm$ideol,ewm$demstr,ewm$legideol), y=ewm$policyli, W=ewm.W, mcmc=100000, beta.mean=c(.44,-.32,.48),beta.var=diag(c(1/69.444,1/204.081,1/69.444)))

#obtain results
summary(ewm.1,quantiles=c(.05,.5,.95))$statistics[1,]
summary(ewm.2,quantiles=c(.05,.5,.95))$statistics[1:2,]
summary(ewm.3,quantiles=c(.05,.5,.95))$statistics[1:3,]
summary(ewm.4,quantiles=c(.05,.5,.95))$statistics[1:2,]
summary(ewm.5,quantiles=c(.05,.5,.95))$statistics[1:3,]

summary(ewm.1,quantiles=c(.05,.5,.95))$quantiles[1,]
summary(ewm.2,quantiles=c(.05,.5,.95))$quantiles[1:2,]
summary(ewm.3,quantiles=c(.05,.5,.95))$quantiles[1:3,]
summary(ewm.4,quantiles=c(.05,.5,.95))$quantiles[1:2,]
summary(ewm.5,quantiles=c(.05,.5,.95))$quantiles[1:3,]
var.terms.1<-sqrt(exp(summary(ewm.1,quantiles=c(.05,.5,.95))$quantiles[2:3,]))
var.terms.2<-sqrt(exp(summary(ewm.2,quantiles=c(.05,.5,.95))$quantiles[3:4,]))
var.terms.3<-sqrt(exp(summary(ewm.3,quantiles=c(.05,.5,.95))$quantiles[4:5,]))
var.terms.4<-sqrt(exp(summary(ewm.4,quantiles=c(.05,.5,.95))$quantiles[3:4,]))
var.terms.5<-sqrt(exp(summary(ewm.5,quantiles=c(.05,.5,.95))$quantiles[4:5,]))
var.terms.1[2,2]/(var.terms.1[1,2]+var.terms.1[2,2])
var.terms.2[2,2]/(var.terms.2[1,2]+var.terms.2[2,2])
var.terms.3[2,2]/(var.terms.3[1,2]+var.terms.3[2,2])
var.terms.4[2,2]/(var.terms.4[1,2]+var.terms.4[2,2])
var.terms.5[2,2]/(var.terms.5[1,2]+var.terms.5[2,2])



####Forest Plot Results####
#Cross-Sectional Results
#EWM RESULTS:
#postscript("ewmProperForest.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("ewmProperForest.pdf",width=5, height=5,  pointsize=10)
lin.mean<-c(.7626,.8415,-1.023,.1301,-.225,.786,1.053,.66,.4426,-.3308,.4778)
lin.lower<-c(.5988,.553,-1.312,-.116,-.4911,.6134,.9576,.565,.2154,-.464,.2452)
lin.upper<-c(.9261,1.13,-.7341,.3764,.04044,.9588,1.148,.7549,.6684,-.1977,.7099)
car.mean<-c(0.6042,0.8197,-0.8483,0.06112,-0.1243,0.6984,1.066,0.6653,0.3693,-0.2951,0.4926)
car.lower<- c(0.4055,0.4876,-1.215,-0.2286,-0.4379,0.5009,0.912,0.5221,0.09221,-0.4931,0.2189)
car.upper<-c(0.8021,1.151,-0.478,0.3501,0.1898,0.8953,1.22,0.8086,0.6458,-0.09462,0.7676)
proper.mean<-c(0.5563,0.6414,-1.0518,0.0802,-0.2100,0.4689,1.0474,0.6292,0.2603,-0.3188,0.4634)
proper.lower<- c(0.4124,0.3907,-1.3423,-0.0916,-0.4457,0.3225,0.9825,0.5614,0.1077,-0.4339,0.2649)
proper.upper<-c(0.7076,0.8912,-0.7619,0.2531,0.0273,0.6123,1.1123,0.6976,0.4164,-0.2085,0.6669)

par(omi=c(.2,1.1,.1,.01))
plot(x=lin.mean,y=c(1:11-.1),xlim=c(-1.3,1.3),ylim=c(0,12), xlab="Coefficient", ylab="",axes=F)
axis(1)
axis(2, at=c(1:11), labels=c('Opinion on party elite','Opinion on Dem ID','Party elite on Dem ID','Opinion on strength','Party elite on strength','Dem ID on strength','Party elite on legis lib','Strength on legis lib','Opinion on policy','Strength on policy','Legis lib on policy'),las=1)
box()
abline(v=0,col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(1:11)-.1,y1=c(1:11)-.1, lty=2)
points(x=car.mean,y=c(1:11+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:11+.1),y1=c(1:11+.1), col='red')
points(x=proper.mean,y=c(1:11+.3), col='blue',pch=3)
segments(x0=proper.lower,x1=proper.upper,y0=c(1:11+.3),y1=c(1:11+.3), col='blue',lty=3)
legend(x=-1.3, y=2.5, pch=c(2,1,3), lty=c(1,2,3), col=c('red', 'black','blue'), legend=c('CAR', 'Independent','Proper CAR'))
#dev.off()
